import math
import numpy as np
import matplotlib.pyplot as plt


plt.rcParams['font.sans-serif'] = ['SimSun']  # matplotlib中文支持，宋体，确保系统字体内有“宋体”
plt.rcParams['axes.unicode_minus'] = False  # 中文字体下的负号支持

def mogi_vertical(x,a,d,g,delta_p,sigma):
    # Calculate the vertical displacement (m)
    # at the surface due to a subsurface
    # pressure source using Mogi's solution
    # x is horizontal distance (m), radial from source
    # a is source radius (m) must be << d
    # d is depth (m)
    # g is shear modulus (MPa), typically 3e5 or less
    # delta_p is the excess pressure of the source (MPa)
    # sigma is Poisson's ratio (typically 0.25)
    w = (1-sigma)*a**3*d*delta_p/(g*(x**2+d**2)**1.5)
    return w  # m displacement

def mogi_horizontal(x,a,d,g,delta_p,sigma):
    # calculate the horizontal displacement (m) at the
    # the surface due to a subsurface pressure
    # source using Mogi's solution.
    # x is horizontal distance (m), radial from source
    # a is source radius (m) must be << d
    # d is depth (m)
    # g is shear modulus (MPa), typically 3e5 or less
    # delta_p is the excess pressure of the source (MPa)
    # sigma is Poisson's ratio (typically 0.25)
    u = (1-sigma)*a**3*x*delta_p/(g*(x**2+d**2)**1.5)
    return u # m displacement
def mogi_V3D(x,a,d,g,delta_p,sigma):
    """计算三维垂直位移并绘图"""
    x=x*10 #单位0.1km
    image_V=np.zeros((2*x,2*x))
    for i in range(2*x):
        for j in range(2*x):
            L_dis=math.sqrt((x-i)*(x-i)+(x-j)*(x-j))*100
            image_V[i,j]=mogi_vertical(L_dis,a,d,g,delta_p,sigma)*1000
    plt.imshow(image_V,origin="lower",cmap="jet_r")
    plt.title("垂直位移图(mm)")
    plt.xlabel("距离(km)")
    plt.ylabel("距离(km)")
    plt.xticks(np.arange(0, 2 * x, int(2 * x / 8)), abs(np.arange(-x/10, x/10, int(2 * x / 80))))
    plt.yticks(np.arange(0, 2 * x, int(2 * x / 8)), abs(np.arange(-x/10, x/10, int(2 * x / 80))))
    plt.colorbar()
    plt.tight_layout()
    plt.show()
def mogi_H3D(x,a,d,g,delta_p,sigma):
    """计算三维水平位移并绘图"""
    x=x*10 #单位0.1km
    image_H=np.zeros((2*x,2*x))
    for i in range(2*x):
        for j in range(2*x):
            L_dis=math.sqrt((x-i)*(x-i)+(x-j)*(x-j))*100
            image_H[i,j]=mogi_horizontal(L_dis,a,d,g,delta_p,sigma)*1000
    plt.imshow(image_H,origin="lower",cmap="jet_r")
    plt.title("水平位移图(mm)")
    plt.xlabel("距离(km)")
    plt.ylabel("距离(km)")
    plt.xticks(np.arange(0, 2 * x, int(2 * x / 8)), abs(np.arange(-x/10, x/10, int(2 * x / 80))))
    plt.yticks(np.arange(0, 2 * x, int(2 * x / 8)), abs(np.arange(-x/10, x/10, int(2 * x / 80))))

    plt.colorbar()
    plt.tight_layout()
    plt.show()





a = 2500.0 #source radius (m)
d = 10000.0 #source depth (m)
sigma = 0.25 #Poisson ratio
delta_p = -400.0 #excess pressure (MPa)
g = 30000.0 #rigidity of elastic half-space (MPa)

# Sakurajima displacements following 1914 eruption
# The data are taken from Mogi (1958)
# The data were collected by line levellin by
# F. Omori and colleagues
# displacement is the difference between 
# measured elevations in 1893 and 1914.
sakura = [(17.0,-291),
(16.2,-291),
(15.75,-302),
(13.5,-407),
(12.45,-446),
(10.4,-527),
(8.7,-658),
(7.85,-770),
(6.7,-894),
(7.75,-776),
(8.25,-703),
(9.15,-608),
(9.85,-526),
(9.7,-536),
(9.55,-560),
(9.3,-577),
(8.5,-684),
(8.55,-690),
(8.95,-681),
(9.15,-685),
(10.05,-616),
(11.1,-481),
(11.9,-427),
(12.15,-369),
(12.1,-348),
(12.5,-301),
(12.3,-306),
(10.95,-433),
(11.05,-426),
(11.0,-453),
(12.0,-356),
(12.4,-344),
(12.5,-363),
(11.6,-401),
(11.2,-468),
(12.6,-379),
(13.95,-314),
(15.8,-242),
(12.15,-369),
(12.5,-338),
(13.75,-260),
(14.15,-244),
(15.4,-206),
(17.25,-155),
(18.8,-143)]




"""
# calculate vertical displacement for this
# range of radial distances (m)
x1 = np.arange(0,40000,1)
x=40
mogi_V3D(x,a,d,g,delta_p,sigma)
mogi_H3D(x,a,d,g,delta_p,sigma)
# plot vertical displacement calculated with
# Mogi solution in mm vs. km
plt.plot(x1/1000, mogi_vertical(x1,a,d,g,delta_p,sigma)*1000)
# Compare the Mogi solution to observations
# for Salurajima following deflation after
# the 1914 eruption
radial, vert_displace = zip(*sakura)
plt.plot(radial, vert_displace, "bo")
plt.xlabel("半径 (km)")
plt.ylabel("垂直位移 (mm)")
plt.title("Mogi模型垂直位移图")
plt.show()
"""